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Abstract 



^ Review of implicit methods of integrating system of stiff ordinary differential 

equations is presented. Defines and graphically presents absolute stability 

<^ region for Gears methods (backward differentiation formula) used to solve 

system of stiff ordinary differential equations. Recommendations for selecting 
(-H the order of Gears method are given. 
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^ 1. Introduction 

^ Ordinary differential equations (ODE) are widely used for modelling real 

processes. Practice shows that the initial value problem (Cauchy problem) 
y—i for ODE systems can be attributed to the following types: soft, stiff, ill- 

2 conditioned and rapidly oscillating. Each type is connected with specific 

T— I demands to integrating methods. Stiff systems can be exemplified by prob- 

j> lems of chemical kinetics [H [2] , nonstationary processes in complex electric 

circuits [3 H], systems emerging while solving equations of heat conduction 
^ and diffusion [S], movement of celestial bodies and satellites [B], plasticity 

^ physics [7j, etc. 

The numerical solution of ODE systems is accompanied by problems due 
to the fact that at the modelling of a complex physical process speeds of local 
processes may vary significantly, while variables systems may be of different 
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orders or change within the interval of integration by orders of magnitude 

Besides, in physical experiments a boundary layer may be observed which 
is characterized by quick changes in the object of research. Such a layer does 
not necessarily emerge at the beginning of the experiment, but only when 
some controlling parameter sharply changes or achieves a certain critical 
value. For the numerical solution of such problems the numerical method 
has to be chosen very carefully [H El I9HT3] . 

The purpose of this paper is to report on a continuing research effort 
aimed at the use and development of numerical methods in computer program 
for plasticity physics. The previous results were reported in [S] while using 
the same mathematical basis for solving stiff ODEs. The paper is organized 
as follows: the review of definitions stiff ODEs presented in Section |2] The 
approaches to finding a numerical solution of the stiff ODEs are discussed in 
Section 3. In Section 4, defines and graphically presents absolute stability 
region for backward differentiation formula and recommendations for various 
types of ODEs are given. 

2. Stiff ordinary differential equations 

Interest for stiff systems appeared at the beginning of the 20th century, 
initially in radio engineering (van der Pol problem, 1920 [51 1^ [TIHTT] ). Then 
there was a new wave of interest in the middle 1950s with application in 
studying equations of chemical kinetics, movement of celestial bodies [6], 
which contained both very slowly and very rapidly changing components. 
Methods like Runge-Kutta, which had been considered as highly reliable, 
produced mistakes in solving such problems. 

One of the first attempts to give a definition of stiff systems was made 
by CF. Curtiss and J.O. Hirschfelder in 1952. They proposed the following 
interpretation: stiff equations are equations where certain implicit methods 
perform better, than using classical explicit ones like Euler or Adams methods 

m- 

There are a number of interpretations of stiffness, each of which reflects 
certain aspects of the numerical solution (e.g. impossibility of using ex- 
plicit methods of integration [19], presence of rapidly damped disturbances 
[m |20] , large Lipschitz constants or logarithmic norms of matrices [211 ES] , 
big difference between eigenvalues of Jacobian matrix [H [21 [231 [2l] , fullness of 
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Jacobian matrix a priori fixed sign of tlie solution number of tran- 
sient pliases [27], etc.) In some applications an important factor influencing 
the behavior of the numerical solution is the order of the ODE system [28j or 
the presence of a boundary level [2H |29] , in others - only the limit behavior 
(gradual change) at wide intervals of integration. It is often unclear whether 
stiffness is attributed to a particular solution or the problem in general. 

In papers P-[3l HSj |23l EH EQ] authors find it difficult to clearly define a 
stiff system ODE because of its complex character, so they present a working 
description of a stiff problem. This is a problem modelling a physical pro- 
cess, components of which possess incommensurable characteristic times, or 
a process characteristic time (reciprocal quantities of Jacobian eigenvalues) 
of which is much smaller than the interval of integration. 

In 1970s L.F. Shampine and CW. Gear, who had gained a wide experi- 
ence of computing experiments with systems having components of the deci- 
sion vector of different orders, offered their own definition of a stiff ODE sys- 
tem: the initial value problem for ODEs is stiff if the Jacobian Jj j = dfi/dz/j, 
i, j = 1, . . . ,N has at least one eigenvalue, for which real part is negative 
with high modulus, while the solution within the major part of the interval 
of integration changes slowly [151 EI] • 

In papers [231 129] ODE system is called stiff if real components of all 
the eigenvalues of Jacobian are negative, i.e. -Re(Aj) < 0, i = 1, . . . , (the 
system is asymptotically stable) and the ratio 

_ max{\Re{\i)\,i = 1,2,...,N} 
^ ~ min{\Re{Xi)\,i = 1,2, ...,iV} 

is large. The parameter s is called the stiffness ratio. 

The problem of defining a stiff system is that for the stiffness ratio the 
boundary value where it becomes big is not given. A system of equations can 
be regarded as stiff if the stiffness ratio s exceeds 10, but in numerous applied 
problems this parameter reaches 10^ and higher [32] . In paper [33] a concept 
of superstiffness is introduced, when the stiffness ratio reaches 10^ . . . 10^^. 

It must be pointed out that there are no simple methods for evaluating 
stiffness, so numerical methods working without stiffness testing are neces- 
sary. 
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3. Approaches to finding a numerical solution of the stiff ODEs 

Let us regard a Cauchy problem for the ODEs of the first order, which 
can be presented as follows: 

Y' = F{x,Y),xe [xo,b],Y{xo) = Fo- (1) 

For numerical integrating of system ([T]) methods using a linear combi- 
nation of the decision vector Yi, Y2, . . . , Yn, . . . and its derivatives in some 
sequence of independent variable Xi, X2, ■ ■ ■ , Xn, ■ ■ ■ are widely applied. Such 
methods are called linear methods [U HI 1^ . 

Since the discovery of the stiffness in the development of numerical meth- 
ods for integrating stiff systems, the following trends have appeared: the 
investigation of stiffness and the establishment of the theoretical apparatus 
stability analysis methods |31 HH], design and enhancement of the methods 
taking into account the specifics of the tasks [SIHIU] and the prospects for 
parallel [HI US] . The most complete overview of the current numerical meth- 
ods for solving stiff ODE systems with an extensive bibliography is presented 
in the papers [3l H Hg [291 ill E] . 

3.1. Taylor series 

One of the easy ways to construct the solution to the system ([T]) at point, 
Xn+i, if it is known at point, method based on the expansion of 

solutions Y{xn+i) in a Taylor series in the neighborhood of point, Xn'- 

y{Xn+l) = yiXn) + hF{Xn, Y^, h), 

where 

F{Xn, Yn, h) = Y'iXn) + hY"{Xn)/2\ + h^Y"'{Xn)/S\ + .... 

If this series is truncated at q-th term and replace Y{xn) with the approximate 
value of Yn, thus an approximate formula is obtained: 



Yn+l = Yn + h{F{Xn, Yn) + hF\Xn, Yn)/2\ + 
+ h^F"{Xn, Yn)/^\ + ... + h'^F^'^^ (Xn, + 1)!). (2) 

If g = 1, the computational scheme of explicit Euler method [H |2l [131 1231 [15] 
is obtained 

Yn+l = Yn + hF{Xn,Yn). 
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Application of the formula is limited to only those tasks which can 
easily calculate the higher-order derivatives of the function F{x, Y) of the 
right side of the system ([T]). Though, it is usually not so. 

3.2. Runge-Kutta methods 

S. Runge (1895), K. Heun (1900) and M. Kutta (1901) put forward an 
approach based on constructing of the formula for Yn+i [H El H]: 

Y{xn+i) = Y{xn) + /i$(x„, F„, h), 

where h - integration step. The function $(■) is close to F{-), but does not 
contain the derivatives from the right side of the equation. Thus, the series 
of explicit and implicit methods requiring s-stage calculation of the right side 
function at each integration step are obtained (s-stage methods). 

The formulas of these methods are ideally applicable for practical calcu- 
lations: they allow to change the integration step h easily. Perhaps the most 
famous is the formula of the 4-th order of 4-stage Runge-Kutta [4] . 

One of the major problems associated with the use of Runge-Kutta meth- 
ods (in fact, almost of all the explicit methods) lies in choosing the size of the 
integration step h, which provides the stability of the computational scheme 
[H El IIHI ESI EH HHl US] . Nevertheless, even nowadays, the explicit adaptive 
methods for solving stiff ODEs [38^ are developed and widely used. 

3. 3. Backward differentiation formulas 

These stiff tasks have made the implicit computational scheme particu- 
larly attractive and led to the development of such implicit methods which do 
not involve calculations based on the size of the integration step [iHl EH [3T| 
H71 - I5T] . The most common among them are the methods of Adams-Moulton 
and "backward differentiation formulas" (more commonly known as Gear 
method). Having got the approximation to the solution at points, Xi, X2, 
. . . , Xn, it is possible to find solutions at point, Xn+i- 

The computational schemes of the Adams-Moulton implicit methods take 
the following form [D El ttSl ESI El SSI 06] : 

Yn-i+l), (3) 

i=0 

where q determines the order of the method, the constants i = 0,1, . . . ,q 
correspond to the chosen order of the method [H EH EI]. Moreover, the 
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implicit Euler's method (the first order) and the trapezoidal method (the 
second order) are the special cases of the last computational scheme ([s]) 
where g = and q = 1- 

The construction of the multistep methods is based on the polynomial 
of the degree q. The approximate value of the solution Y{x) at the point, 
Xn+i appears as a linear combination of the several approximate values of 
the solution and its derivative at this and the previous q points. Obviously, 
the use of the multistep formulas requires the calculations of the q units of 
the initial values Yi, I2, • • • , Yn- The accuracy of setting these q units should 
not be less precise than the accuracy of the formula. Thus, the polynomial 
can be represented PQ [21 HSl [221 [211 [ISl [IS] as the following formula: 

g g 

Yn^j + l). (4) 

i=l j=0 

Some constants and (3i in equation ^ can take zero values. When (3i = 
= /32 = ■ ■ ■ = Pq = 0, it is possible to construct backward differentiation 
formulas. 



4. The region of implicit methods stability 

The first studies on the stability of multistep methods refer to the re- 
searches of G. Dahlquist [521 [53]. According to the definitions above, the 
stiff ODEs require the stability of numerical methods used to solve them. 
When getting the asymptotically stable solution of the stiff Cauchy problem, 
an error of the difference method should non-increase under any step, i.e. the 
method should be absolutely stable. Current reviews of the stability regions 
of multistep methods can be found in |3l [1] . 

The researchers [231 [13 [51] give a more clear definition to the stability 
method through the model first-order equation 

y' = >^y,yixo) = yo- (5) 

The general solution of the equation ([s]) takes the form of 

y = C ■ exp{\x}, 

where C - a constant, the solution of the equation ([s]) - a function y = 
yo ■ exp{\{x — Xq)} - tends to zero, if Re{X) < and it infinitely grows in 
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absolute value for Re{X) > 0, where A is a complex number. Further, the 
concept of absolute stability and the study of the absolute stability of the 
numerical methods are considered for the model equation dsl). 

Here, the stability domain of the multistep method ffi of solving the 
initial value problem ([s]) is defined as the set of points of complex numbers 
plane defined by the complex variable a = hX. For a = hX, this method 
applied to the model equation ([s]) is stable, i.e. it provides non-increase of 
an error [23l SSI IM] • 

To determine the region of the implicit methods stability Q, the charac- 
teristic (complex) polynomial is used: 

P{Z) = {Z'i - aiZ'i-^ - 02-2""^ - ... - aq) + 

+a{Poz'' + Piz'i-^ + /32^^-2 + . . . + = 0, 
it may be represented in the form of [3l HI [6l [211 EI] '■ 

a{e) = - '=° , (6) 

where a^, (3^ - coefficients of the method {k=Q, 1, 2,. . . , q; f3i = = . . . = 
= /3q = 0), g - the order of the method, z - a unit imaginary number, z = e*^ 
- an imaginary number, < 6 < 2tt. To determine the region of the method 
absolute stability for the given value cr = ctq, it is necessary to find the 
solution of the equation ^ relatively to 9. 

A set of points generated by the equation ^ corresponds to a geometrical 
locus of points of single radicals T^, for which \z\ = 1 is true. The region of 
absolute stability of implicit methods Q is considered to be an external area 
Fo- since at \a\ = h\X\ — )■ oo imphcit methods Q are stable [3t HI [2H l3T]. 

To determine a geometrical locus of points described in the equation (6), a 
computer algebra system MathCAD is used. As a result, a geometrical locus 
of single radicals F^. is obtained (fig. [T]). The region of absolute stability of 
methods is considered to be an external region F^. (shaded). 

The study region of stability (F^- is a simple closed curve, fig. [T]) shows 
that implicit methods ^ from 1 to 6 order inclusive are stiffly stable (first 
introduced by [31] )• The stiffness property of the implicit method ^ from 
1 to 6 order inclusive is attained at various values of a real number 6 < 
(fig. [T| c — (7 marked with a dotted line). Particularly, for the methods of 
the first and second order it is 5 = —0.; for the third order 6 = —0.1; for 
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the fourth order 6 = —0.7; for the fifth order 5 = —2.4, for the sixth order 
6 = -6.1. 

For Gear's method of the seventh order, the equation (|6| will be: 



ai9) 



m _|_ 980 j66» _ 490 ise _ 4900 i4e_ 

~'~ 363 121'^ 1089 

1225 i36> _|_ im 126 _ 490 ie , 20 

363 "'"121^ 1089 "'" 363 

140pi7e 
369^ 



Finding a solution to a{d) relative to 9, it turns out that Gear method of 
the seventh order does not meet the requirements of stiff stability (fig. 2)|24[ 
121]. At origin of coordinates and at Re{X) ~ —8, there are intersection points 

of r,. 

The results obtained for the equation ([s]) can be distributed on the ODEs 
|44j . In the case of an autonomous system of ODE of Y' = AY type, where 
A is a constant matrix, it becomes possible to transact matching Jordan's 
form and proceed to look for a solution to the system of ODE of Z' = 3Z 
type, where J = T^^AT = diag{Xi, A2, . . . , A„}, Aj - eigenvalues of matrix 
A, i = 1,2,...,A^, r = TZ, Z = T'^Y. The matrix T is composed of 
eigenvectors of the matrix A. Thus the initial system of ODE decomposes 
into n scalar equations, for which the solution can be found and the above 
approach to the region of stability determination applied [H [211 EI] . 

If the coefficients of the system Y' = A{x)Y are not constant, the check 
of the eigenvalues A at each value x becomes laborious to calculate [31]. It 
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should be noted that operating with nonhnear systems of Y' = AY + 
+ G{x, y) type, the stabihty of the solution can be provided only at the origin 
of coordinates, moreover the stability can be broken for eigenvalues located 
on an imaginary axis |T9]. 

The overview of the alternative ways of defining stability regions for the 
implicit methods can be seen in HI [12] . 

Implicit methods (Gear methods) can be applied for the calculation of 
a big category of the stiff ODEs. In this case decrease of stepsize (to the 
minimum possible) doesn't always let us adapt to local solution and decrease 
volume of computation with required precision. Optimal strategy of using 
multistep methods implies availability of order autocontrol (from 1 to 6) and 
stepsize. 

5. Conclusion 

The area of the numerical methods for solving of ODEs is one of the 
most well-investigated topics in the mathematical literature. A number of 
techniques and solvers have been suggested, development, and described, but 
a clear definition of stiffness has not been provided, so the working definition 
of stiffness is still topical. 

As a rule in most solvers for ordinary differential equations the explicit 
first order method (Euler method) or a second order method (trapezoidal 
method) is applied. Implicit Gear methods (backward differentiation for- 
mulas) are stiff from 1 to 6 order inclusive, so for the acceleration of the 
integration process of ordinary differential equations increasing order could 
be applied. 

The results of the calculations let us define absolute stability regions for 
the implicit methods where changing of integration step over wide region 
when computational stability of the method is constant. 
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